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1 Introduction 

Formal Laurent-Puiseux series of the form 

CO 

f{x) = £ a k x k / n (1) 

k=ko 

with coefficients a k £ C (fc € 2) are important in many branches of mathematics. Maple supports 
the computation of truncated series with its series command, and through the powerseries 
package || infinite series are available. In the latter case, the series is represented as a table of 
coefficients that have already been determined together with a function for computing additional 
coefficients. This is known as lazy evaluation. But these tools fail, if one is interested in an explicit 
formula for the coefficients ak- 



in this article we will describe the Maple implementation of an algorithm presented in ||-|13] 
which computes an exact formal power series (FPS) of a given function. This procedure will enable 
the user to reproduce most of the results of the extensive bibliography on series @. We will give 
an overview of the algorithm and then present some parts of it in more detail. 

This package is available through the MAPLE-share library with the name FPS. We will flavor 
this procedure with the following example. 
> FormalPowerSeries (sin(x) , x=0) ; 

infinity 

k (2 k + 1) 

\ (-1) x 
) 

/ (2 k + 1) ! 

k = 



2 Preliminary Results 



To deal with many special functions, it is a good idea to consider the (generalized) hypergeometric 
series 



b\ b 2 



On 



_ \ - (ai)fc • (ag)fc • • • (aph k 
t (bi)k-(b 2 ) k ---(b q ) k k\ X 



(2) 
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where {a)k denotes the Pochhammer symbol (or shifted factorial) defined by 



(o)jfc 




if k = 
if fc G N 



We have defined this function in our package under the name pochhammer. 

oo 

The coefficients of the hyper geometric series J2 A^x are the unique solution of the special 

k=0 

recurrence equation (RE) 

a - (fc + gi) • (k + a 2 ) • • • (k + op) 

^ fc+1 : " (k + 60 • (fc + 6„) • • • (Jb + b q )(k + l)' Ak {ke ^ } 

with the initial condition 

A :=l. 

Note that ^ +1 is rational in fc. Moreover if ^| +x is a rational function i?(/c) in the variable k then 
the corresponding function / is connected with a hypergeometric series; i.e., if k = —1 is a pole 
of R, then / corresponds to a hypergeometric series evaluated at some point Ax (where A is the 
quotient of the leading coefficients of the numerator and the denominator of R); whereas, if k = — 1 
is no pole of R, then / may be furthermore shifted by some factor x s (s £ Z). 

We further mention that the function / corresponding to the hypergeometric series 



fix) := p F t 



a\ a 2 • • • a p 
6i b 2 ■■■ b q 



satisfies the differential equation (DE) 

0(0 + b x - 1) • • • (9 + b q - l)f = x(9 + 01 ) • • • (9 + op)/ (3) 

where 9 is the differential operator x-^ . An inspection of the hypergeometric DE (|3|) shows that 
it is of the form 

EEV/ 0) = o (4) 

3=0 1=0 

with certain constants € C and Q = max{p, q} + 1. Because of their importance in our devel- 
opment, we call a DE of the form (|j), i.e. a homogeneous linear DE with polynomial coefficients, 
simple. 

We extend the considerations to formal Laurent- Puis eux series (LPS) with a representation 

oo 

/:= E a ^ k/n K + 0) (5) 

k=ko 

for some ko € 2, and n 6 IN. LPS are formal Laurent series, evaluated at yfx. A formal Laurent 
series (n = 1) is a shifted FPS, and corresponds to a meromorphic / with a pole of order — k$ at the 
origin. The number n in development (]|) is called the Puiseux number of (the given representation 
of) /• 
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Definition 1 (Functions of hypergeometric type). An LPS f with representation ^) is 
called to be of hypergeometric type if its coefficients a k satisfy a RE of the form 

a k+m = R(k) a k for k > k (6) 
a k = A k for k = k ,k + 1, . . . , k + m - 1 

for some m E N, A k £ (D (k = ko + 1, k$ + 2, . . . , ko + m — 1), A ko 6 C \ {0}, and some rational 
function R. The number m is then called the symmetry number of (the given representation of) f . 
A RE of type (^) is also called to be of hypergeometric type. 

We want to emphasize that the above terminology of functions of hypergeometric type is pretty 
more general than the terminology of a generalized hypergeometric function. It covers e.g. the 
function sinx which is not a generalized hypergeometric function as obviously no RE of type (^) 
with m = 1 holds for its series coefficients. So sinx is not of hypergeometric type with symmetry 
number 1; it is, however, of hypergeometric type with symmetry number 2. A more difficult example 
of the same kind is the function e aTCSmx which is neither even nor odd, and nevertheless turns out to 
be of hypergeometric type with symmetry number 2, too. Also functions like x~ 5 since are covered 
by the given approach. Moreover the terminology covers composite functions like sin ^/x, which do 
not have a Laurent, but a Puiseux series development. Each LPS with symmetry number m, and 
Puiseux number n, can be represented as the sum of nm shifted m-fold symmetric functions. 

We remark further that one can extend the definition of functions of hypergeometric type to 
include also the functions of the form J f where / is an arbitrary LPS (see ||, Section 8). Note 
that because of the logarithmic terms these functions, in general, do not represent LPS. 

In [|| it has been proven, that each LPS of hypergeometric type satisfies a simple DE, which 
is essential for our development. Now assume, a function / representing an LPS is given. In order 
to find the coefficient formula, it is a reasonable approach to search for its DE, to transfer this 
DE into its equivalent RE, and you are done by an adaption of the coefficient formula for the 
hypergeometric function corresponding to transformations on / which preserve its hypergeometric 
type. Below, we will discuss these steps in detail. 

The outline of the algorithm to produce a formal Laurent-Puiseux series expansion of the 
function / with respect to the variable x around x = is shown in Algorithm [I]. Note that series 
expansions around other points can easily be reduced to this case. 

3 Search of the DE 

In this section we present the algorithm that searches for a simple DE of degree k for a given 
function /. We set up the equation 

fW( x ) + J2A j fW(x) = 

3=0 

and expand it. Then we collect the coefficients of all the rationally dependent terms and equate 
them to zero. For testing whether two terms are rationally dependent, we divide one by the other 
and test whether the quotient is a rational function in x or not. This is an easy and fast approach. 
Of course, we could also use the Risch normalization procedure to generate the rationally 

independent terms, but first this normalization is rather expensive and only works for elementary 
functions and second, our simplified approach will not lead to wrong results, it may at most happen 
that we miss a simpler solution, which, in practice, rarely happens, however. This procedure by its 
own is available under the name SimpleDE. 



3 



FormalPowerSeries (f , x=0) 
for k := 1 to k max do 

Search a simple DE of degree k of the form 

DE : = f( k )(x) + J2A j fW(x)=0 (7) 

where the Aj are rational functions in x. 
if the search was successful, then 

Convert the DE into a recurrence equation of the form 

M 

RE := ^Pja k+j = (8) 

3=0 

for the coefficients where pj are polynomials in k and M £ IN 
if the RE only contains one or two summands then 

/ is of hypergeometric type and the RE can be solved 
elif the DE has constant coefficients then 

/ is of exp-like type and the RE can also be solved. 

fi 

fi; 

if f^ k \x) is a rational function in x, then 

use the rational Algorithm and integrate the result k times. 

f i 



od 



Algorithm 1: Algorithm FormalPowerSeries 



The resulting differential equation only depends on the form of the derivatives f^\x) (which 
of course must be known to the system). As an example we look at the Airy wave function Ai, 
whose derivative presently (in Maple) is given as 

> diff (Ai(x) ,x) ; 

1/2 3/2 
x 3 BesselK(2/3, 2/3 x ) 

- 1/3 

Pi 

The simple DE then becomes 

> SimpleDE(Ai(x) ,x) ; 

/ 3 \ / 2 \ 

Id I 2 / d \ I d I 

I F(x) I x - x I F(x) I - I F(x) I = 

13 I \ dx / I 2 I 

\ dx / \ dx / 

which is not the simple DE of lowest degree valid for Ai. This happens since the second derivative 
is not expressed in terms of Ai (and diff (Ai (x) ,x)) itself (or in other words since Ai itself is not 
expressed in terms of BesselK). We may introduce a new function newAi by defining its derivatives 
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> 'diff/newAi 1 := (e,x) -> diff (e,x)*newAiPrime(e) : 

> 'dif f /newAiPrime' := (e,x) -> dif f (e ,x) *e*newAi (e) : 

> SimpleDE(newAi (x) ,x) ; 

/ 2 \ 
I d I 

I F(x) | - x F(x) = 

I 2 I 
\ dx / 

and we get the expected differential equation for the Airy wave function Ai. 

It may happen, that for a given function f(x) a DE of degree k exists, but which has neither 
constant coefficients (which we call the explike case) nor is the corresponding RE of hypergeometric 
type and hence no closed form for the LPS can be computed. What we can do in this situation is 
to look for a DE of higher degree which then will have free parameters as from the existence of a 
DE of degree k follows the existence of families of DEs of higher degree. These parameters can be 
set freely and we can try to choose them in such a way, that we can use our tools to compute a 
formal LPS, i.e. we either need a RE of hypergeometric type or a DE with constant coefficients. 

For the first case we convert the DE into the corresponding RE of the form @ and try to set 
all the coefficients but two to zero. For example, let f(x) = x _1 e x sinx. We find DEs of degree 2 
and of degree 3, but none of the corresponding REs can be solved. The RE which corresponds to 
the DE of degree 4 has the form 

4 

^2pj{k)a k+j = 

3=0 

where 

p (k) = 2A 2 +4A 3 + A 

Pl (k) = -AA 2 - 10A 3 - 16 - 2A 3 k - 2A 2 k 

p 2 {k) = 18A 3 + 5A 2 k + 6A 2 + 48 + 8k + 6A 3 k + A 2 k 2 

p 3 {k) = (4 + k)(A 3 k 2 + 2A 3 k - 24 - 3A 3 ) 

Pi {k) = (k + 5)(4 + k)(k 2 + k + 6) 



The solution of the equations pi(k) = 0,p 2 (k) = 0,p 3 (k) 
sum remain) with respect to A 2 and A 3 is 

„ k 2 + 5k + 12 



A 2 



k 3 + 4k 2 + k-6' 



A, 



(forcing that only two terms of the 
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k 2 + 2k-3' 



This results in the following RE of hypergeometric type (after multiplying by ^jjp^fjp' ) 

(k + 5)(k + 4) (k + 3) (k + 2) a k+1 + Aa k = 

with symmetry number m = 4. Of course, we try to keep the symmetry number, that is the number 
of resulting sums, as small as possible. The final result for this example is 
> FormalPowerSeries(exp(x)*sin(x)/x, x) ; 



/ infinity 



\ 

) 

/ 

\ k = 



k (- 
(-1) 64 



\ 

k) k (4 k) I 
256 x I 



(1 + 4 k) ! 



/infinity 

■ k ( 

\ (-1) 64 
) 

/ (1 + 4 k) ! (2 k + 1) 



\ 

k) k (1 + 4 k) I 
256 x I 



/ \ k = 



5 



/infinity \ 

' k (- k) k (4 k + 2) I 

\ (-1) 64 256 x I 
) 1 

/ (4 k + 3) (1 + 4 k) ! (2 k + 1) I 



\ k = / 

Note that one of the four sums (with the powers x 4k+3 ) vanishes as a result of a vanishing initial 
coefficient. 

If this step fails, then we try to choose the parameters such that the DE gets constant coefficients. 
Let us look at a rather similar example, f(x) = x e x sin(2x). We again find DEs of second and third 
order, but their corresponding REs cannot be solved. The DE of degree 4 has two free parameters 
A 2 and A3 as expected 

ax* ax 6 ax z 

((A 3 - 2 A 2 + 12) x 2 + (-6 A 3 - 2 A 2 + 4) x) -^f(x) + 

((10 A 3 + 5 A 2 - 5)x 2 + (14 A 3 + 2 A 2 + 28)x + (6 A 3 + 2 A 2 - 4)) f(x) = . 

The DE will have constant coefficients, if both A 2 and A 3 are constants and if they meet the 
following equations 

6A 3 + 2A 2 -4 = 
UA 3 + 2A 2 + 28 = ' 

The solution of this system of equations is A 2 = 14 and A 3 = —4. If we insert these values in the 
differential equation and divide by x 2 , then we get the DE 

~j~if( x ) ~ 4^/(x) + U-^f(x) ~ f{x) + 25/(x) = 
ax* dx^ ax* ax 

which has constant coefficients leading to a constant coefficient RE for b). given by f(x) = J2 tt x k 
that can be solved: 

> FormalPowerSeries(x*exp(x)*sin(2*x) ,x) ; 

infinity 

/ k 1/2 k 1/2 \ 

\ I (5 ) cos(k arctan(2)) k (5 ) sin(k arctan(2)) k| k 

) 1-2/5 + 1/5 1 x 

/ \ k! k! / 

k = 



4 Conversion to the Recurrence Equation 

As it has been proven in this transformation is done by the substitution 

x lfi3)( x )^(k+l-l) r a k+j ^ (9) 

into the DE. 

We give here a small Maple procedure which performs this transformation. We assume, that 
the j th derivative of f(x) is represented by the expression f (j). You may compare the procedure 
ConvertDEtoRE with the rule based solution presented in Q. 
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> ConvertDEtoRE := proc(de, f, x, a, k) local X, F, 1, j; 



> if type(de,'+') then 

> map (ConvertDEtoRE, de, f, x, a, k) 

> else 

> X := select (has, j*de, x) ; 

> F := select(has, j*de, f); 

> 1 := degree(X, x) ; 

> j := op(l,F); 

> de/X/F * pochhammer(k+l-l, j) * a(k+j-l) 

> fi 



> end: 

The following example converts the left hand side of the DE for e x into the corresponding left hand 
side of the RE for the coefficients a k of the FPS of e x . 

> ConvertDEtoRE(F(l)-F(0) , F, x, a, k) ; 

(k + 1) a(k + 1) - a(k) 

The search of a DE and its conversion to the RE is directly available through the command 
SimpleRE. We see that newAi(x 2 ) is of hyper geometric type with symmetry number m = 6, whereas 
the second RE is not of hyper geometric type. 

> SimpleRE (newAi(x~2) , x) ; 

(k - 1) (k + 1) a(k + 1) - 4 a(k - 5) = 

> SimpleRE(x/(l-x-x~2) , x) ; 

(1 - k) a(k) + (k - 1) a(k - 1) + (k - 1) a(k - 2) = 

The latter example is the generating function of the Fibonacci numbers, and we get the expected 
RE. Note, that the common factor (k — 1) ensures, that the RE holds V k £ 7L. 

5 Solving a Recurrence Equation of Hypergeometric Type 

If the recurrence equation of a function f(x) is of the form 

Q(k) = P{k)a k (10) 

(P, Q polynomials) then / is of hypergeometric type and the corresponding series representation has 
symmetry number m. The explicit formula for the coefficients can be found by the hypergeometric 
coefficient formula @ and some initial conditions. Based on the analysis of the polynomials P(k) 
and Q(k), we will convert the RE into one corresponding to a Taylor series by applying a sequence 
of transformations stated in the following lemma which preserve the hypergeometric type. 

Lemma 1 Let f be a formal Laurent series (FLS) of hypergeometric type with representation (^), 
whose coefficients a k satisfy a RE of the form ^j, then the following functions are of hypergeometric 
type, too. Their coefficients b k satisfy a RE whose relation to the RE of f is also given. 



(a) 


x n f{x) 


bk+m ~ 


- R{k — n) bk 




(b) 


f{Ax) 


bk+m = 


-- A m R{k) b k 


A e C 


(c) 


f(x n ) 


bk+n m 


= R{k/n) b k 


n G IN 


(d) 


f(x 1/m ) 


bk+1 = 


R(k m) b k 




(e) 


/'(*) 


bk+m = 


--^^R(k + l)b k 
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For a proof of this lemma we refer to ||, Lemma 2.1 and Theorem 8.1. 

First of all, we inspect the roots of P(k) and Q(k) of the RE (|To|). If there are any rational 
roots, then we know that / corresponds (possibly) to a Puiseux series. In this case we transform / 
to a function of Laurent type by an application of transformation (c) , where n is the least common 
multiple of the denominators of all rational roots of P(k) and Q(k). The FLS we get when we solve 
the transformed RE can be transformed back to the LPS of / by substituting x by x 1//n . 

Let us now assume that / can be expanded in a FLS. We reduce this problem to solving the RE 
of a function which has a FPS expansion. For that we remove the finite pole of / at the origin by 
multiplying / with a suitable power of x (transformation (a)). From the information of the RE we 
can determine which power we have to use. Let us assume that / may be expanded in a Laurent 
series, i.e. that 3 k$ : V k < k^ : a k = 0. From these known coefficients we can derive further 
ones using the given RE in the form 

P(k) . . 

ak+m= Q(k) ak - 

If we know that a k = then also ak+ m = provided that Q{k) ^ 0. Let k m i n be the smallest 
integer root of Q{k). Consequently at = V k < k m i n + m. From this it follows, that 

g — x —{k m in+m) J 

may be expanded into a FPS, i.e. g has no pole at the origin, given the assumption that / may be 
expanded in a Laurent series. This latter assumption can be tested by computing the limit of g as 
x tends to 0. This limit must be finite. (Since we also allow logarithmic singularities, we test in fact 
whether the limit of x ■ g is 0.) If this is not the case, then the assumption that / may be expanded 
into a FLS is wrong and / must have an essential singularity, i.e. ^ ko : V k < ko : a k = 0. 
Otherwise the FPS of g exists and can again easily be transformed back to the FLS of /. 

What we finally have to show is how to solve a RE which corresponds to a given function / 
which has a FPS expansion. The RE (|ll|) is valid V k : Q(k) ^ 0, especially V k > k max where 
kmax is the largest root of Q(k). Consequently we must determine a k for k = 0, 1, ... , k max + m 
and have to solve the hypergeometric RE for x~ k /(x 1//m ) using (Q) for k > k max . We investigate 
now how this condition can be weakened. 

The coefficient a k is given by the limit lim f^ k \x)/k\. Let us assume, that this limit is finite. 

a:— >0 

Then the hypergeometric RE can be solved in the case that V j > : Q(k + j m) ^ 0, otherwise 
simply x k is added to the result. Moreover note, that 

(P(k) = V a k = 0) A V j > : Q(k + j m) ^ A a k is finite =^ V j > : a(k + j m) = 

and in this case the RE does not have to be solved and it is enough to add a k x k to the result. The 
indices k + j m, j > no longer need to be considered. 

If dfc is infinite, then we found a logarithmic singularity which we can remove by working with 

g=(x- k f(x))> 

and by integrating and shifting the resulting power series S. The RE of g can be obtained from 
the RE of / by applying transformations (a) and (e). The constant term which we lose by the 
differentiation can be determined by computing the limit 

lunf(x)/x k - [S(x)dx. 
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Note that g in general is a function with a FLS expansion and we first must remove the pole to get a 
function with a FPS expansion. This is the reason why we have chosen a recursive implementation 
of the RE solver. The procedure hypergeomRsolve accepts as parameters f(x),P(k),Q(k) and 
the symmetry number m. Every application of a transformation of Lemma |l] is nothing else 
but a recursive call of the RE solver. One may inspect which steps the algorithm performs by 
assigning the variable inf olevel [FormalPowerSeries] . As an example we trace the computation 
for f(x) = x _1 sin yfx. A DE of degree 2 is found whose corresponding RE is of hypergeometric 
type. From the root of Ik + 3 it follows that the Puiseux number is 2 and so transformation (c) 
with n = 2 is applied. The resulting function is of Laurent type. The smallest integer root of Q(k) 
of the transformed RE is —4 and the symmetry number is m = 2, hence we multiply the function 
with x 2 and adjust the RE accordingly. We end up with the function sin a; whose FPS can be 
computed directly. This result is then transformed back to the LPS of f{x) according to the two 
transformations we applied. 

> inf olevel [FormalPowerSeries] := 4: 

> FormalPowerSeries(sin(sqrt(x))/x,x) ; 



FPS/FPS 
FPS/FPS 
FPS/FPS 
FPS/FPS 



looking for DE of degree 1 
looking for DE of degree 2 
DE of degree 2 found. 
DE = 

2 

4 x F" (x) + 10 x F 



FPS/hypergeomRE : 
FPS/hypergeomRE : 
FPS/hypergeomRE : 
FPS/hypergeomRE : 
FPS/hypergeomRE : 

FPS/hypergeomRE : 
FPS/hypergeomRE : 
FPS/hypergeomRE : 
FPS/hypergeomRE : 
FPS/hypergeomRE : 

FPS/hypergeomRE : 
FPS/hypergeomRE : 
FPS/hypergeomRE : 
FPS/hypergeomRE : 
FPS/hypergeomRE : 
FPS/hypergeomRE : 
FPS/hypergeomRE : 



(x) + (2 + x) F(x) = 



RE is of hypergeometric type. 
Symmetry number m = 1 

RE: 2 (k + 2) (2 k + 3) a(k + 1) = - a(k) 
RE modified by k = l/2*k 
=> f := sin(xj/x~2 

RE is of hypergeometric type. 

Symmetry number m = 2 

RE: (k + 4) (k + 3) a(k + 2) = - a(k) 

working with x~2*f 

=> f := sin(x) 

RE is of hypergeometric type. 
Symmetry number m = 2 

RE: (k + 2) (k + 1) a(k + 2) = - a(k) 
RE valid for all k >= 
a(0) = 

a(2*j) = for all j>0. 
a(l) = 1 



infinity 
\ 



k (k - 1/2) 
(-1) x 

(2 k + 1) ! 



k = 



> inf olevel [FormalPowerSeries] 



6 Rational Algorithm 

If the given function (or any of its derivatives) is rational in x we can apply the rational algorithm 
as described in Section 4 of ||. First the complex partial fraction decomposition of f(x) has to be 
calculated. Each term of the form , _f y can be expanded by the binomial series whose coefficients 
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arc 



> FormalPowerSeries(l/( (x-1) ~2*(x-2) ) ,x) ; 
infinity 



l) j c ( j + k-1 
k 



k k k 

\ (k! - 2 2 k! + 2 (1 + k) ! 2 ) x 

) - 1/2 

/ k 

2 k! 

k = 

> FormalPowerSeries((l+x+x~2+x~3)/((x-l)*(x-2)) ,x) ; 

/infinity \ 
| | 

I \ (8 2 - 15) x I 

x + 4 + I ) 1/2 1 

I / k I 

I 2 I 

\ k = / 

> FormalPowerSeries(C/(B*A - A*x - B*x+x~2),x); 

infinity 
jj k k 

\ C(AA-BB)x 
) 

/ k k 

(A - B) B B A A 

k = 

To get the complex partial fraction decomposition we must factor the denominator which may 
be rather complicated, hence the rational algorithm to compute the full partial fraction expansion 
presented in [Q] may be used. The following example uses this code. This method can be forced to 
be used, if the environment variable _EnvExplicit is set to false. 

> FormalPowerSeries(l/(x~4+x+l) ,x) ; 

infinity 

/ 2 3\ 

\ I \ 27+64 alpha - 48 alpha + 36 alpha I k 

) | ) 1/22 9 1 x 

/ I / (1 + k) I 

I alpha I 

k = \alpha = °/„l / 



7.1 : = 



RootOf(_Z + _Z + 1) 



A closed form of the Fibonacci numbers can be derived by computing the FPS of their generating 
function x/(l — x — x 2 ). 

> expand (FormalPowerSeries (x/ (l-x-x~2) ,x)) ; 



/ infinity 



1/5 5 



1/2 



) 



/ k 

I x 

| 

I 1/2 k 

\ (1/2 5 - 1/2) 



\ 

k \| 
x I I 
1 | 

1/2 k| I 

(- 1/2 - 1/2 5 ) /I 
/ 



\ k = 

The expand command converts the coefficient in the usual notation. If the factorization of the 
numerator is avoided, then the following result is obtained: 
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> _EnvExplicit := false: 

> FormalPowerSeries (x/ (l-x-x~2) ,x) ; 

infinity 

/ \ 

\ I \ alpha - 2 

) | ) _ 1/5 

/I / (1 + k) 



k 

x 



I alpha 

k = I 2 

\alpha = RootOf (- 1 + _Z + _Z ) / 

> _EnvExplicit := ' _EnvExplicit ' : 

7 Special Functions, in Particular Orthogonal Polynomials 

The algorithm has been extended to handle many special functions, in particular orthogonal poly- 
nomials [Oj. Our implementation covers this approach. 

We have seen, that the only precondition which must be met by a function to be covered by 
the algorithm is that its derivative is defined in terms of functions which also may be handled by 
our algorithm. 

In the case of families of orthogonal polynomials we have the following special situation: The 
general derivative of an orthogonal polynomial of degree n can be defined in terms of the orthogonal 
polynomials of degree n, and n — 1. But on the other hand, furthermore a recurrence equation is 
known which allows to express the polynomial of degree n in terms of the polynomials of degree 
n — 1, and n — 2. Combining these two facts, it is possible to find a second order DE for the general 
polynomial of degree n. If the resulting RE is of hyper geometric type (which depends on the point 
of expansion) the algorithm further yields an LPS representation for the general polynomial of 
degree n. 

Similarly if a function family F(n, x) possesses a differentiation rule of the form 

^ — - = po(n, x) F(n, x) + pi(n, x) Fin — 1, x) 

ox 

with rational expressions po, and p\ (or a similar rule with m rather than 2 expressions on the 
right) then by the product and chain rules of differentiation the second derivative has the form 

7T ' X ^ = qo(n, x) F(n, x) + qi(n, x) F(n - 1, x) + qi{n, x) Fin - 2, x) 
ox z 

with rational expressions qo,qi, and q2, and all higher derivatives of F(n,x) obtain similar repre- 
sentations. Each differentiation increases the number of representing expressions by one. However, 
if we further know a recurrence equation of the form 

F(n, x) = r\(n, x) F(n — 1, x) + r2(n, x) F(n — 2, x) 

with rational expressions r\, and r<i (or a similar equation with m rather than 2 expressions on 
the right) then a recursive application of this equation can be used to simplify each combination of 
derivatives of F(n, x) to a sum of 2 (or m, respectively) rationally independent ones. 

To give an example, we consider the Fibonacci polynomials. The recurrence equation of the 
family of Fibonacci polynomials F n (x) is 

F n (x) = xF n _i{x) + F n - 2 {x) 
F (x) = 
Fi(x) = 1. 
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We can teach our procedure to use this recurrence equation by assigning the table FPSRecursion. 
The second index specifies the number of arguments of the function family. 

> FPSRecursion [Fibonacci , 2] := (n,x) -> x*Fibonacci (n-1 ,x) + Fibonacci (n-2 ,x) : 
The derivative rule is given by 

dF n {x) = (n-l)xF n {x) + 2nF n _ 1 {x) 
dx x 2 + 4 

and written in Maple 

> 'diff /Fibonacci' := proc(n, e, x) 

> diff (e , x) * ( (n-1) *e*Fibonacci (n, e) +2*n*Fibonacci (n-1 , e) ) / (e~2+4) 

> end: 

This rule has been derived as follows. First, an explicit formula of F n (x) has been computed 

oo 

using our algorithm to generate the FPS of the generating function F n (x) t n = 1 _ x t t _ t -2 of the 

n=0 

Fibonacci polynomials that is an easy consequence of the recurrence equation: 

> FormalPowerSeries(t/(l-x*t-t~2) , t) ; 

2 1/2 k 2 1/2 k k 

(- (- 1/2 x - 1/2 (x + 4) ) + (- 1/2 x + 1/2 (x +4) ) ) t 

Sum(- , 

2 1/2 k 2 1/2 k 2 1/2 

(- 1/2 x - 1/2 (x + 4) ) (- 1/2 x + 1/2 (x + 4) ) (x + 4) 

k = . . infinity) 

After some simplifications, the coefficient of this FPS, i.e. F n {x) has the following form: 

> F : = ((l/2*x+l/2*(x-2+4)-(l/2))-n-(l/2*x-l/2*(x-2+4)-(l/2))-n)/(x-2+4)-(l/2) ; 

2 1/2 n 2 1/2 n 

(1/2 x + 1/2 (x + 4) ) - (1/2 x - 1/2 (x + 4) ) 



2 1/2 
(x + 4) 

We now make the following ansatz for the derivative rule, namely 

F n (x)' = aF n (x) + bF n _ 1 (x) 

and try to solve this equation for the unknown parameters a, and b: 

> eq := diff(F,x) - (a*F + b*subs (n=n-l ,F) ) : 

> eq := numer (normal (eq, expanded)): 

> indets (eq) ; 

2 1/2 2 1/2 n 

{n, a, b, x, (x + 4) , (1/2 x + 1/2 (x + 4) ) , 

2 1/2 n 2 3/2 

(1/2 x - 1/2 (x + 4) ) , (x + 4) } 

> solve ({coeffs(eq, {"[5.. 8]})}, {a,b}) ; 

x (n - 1) n 

{a = , b = 2 } 

2 2 
x + 4 x + 4 

> assign(") ; 

> a*Fibonacci(n,x)+b*Fibonacci(n-l,x) ; 
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x (n - 1) Fibonacci (n, x) n Fibonacci (n - 1, x) 
+ 2 

2 2 
x + 4 x + 4 

Note that, again, the above procedure essentially equates the coefficients of the rationally indepen- 
dent terms of our setting to zero. 

The algorithm trying to find a DE computes the first and the second derivative of F n (x), 
expresses all occurrences of F n (x) in terms of F n -i(x) and F ra _2(x), and finally returns the following 
solution: 

> SimpleDE (Fibonacci (n, x) ,x) ; 

/ 2 \ 

2 I d I / d \ 

(x + 4) I F(x)| - (n - 1) (n + 1) F(x) + 3 x I F(x) I = 

12 1 \ dx / 

\ dx / 

If we declare the initial value for x = which is for even n and 1 for odd n (which follows from 
the recurrence equation for x = 0: F n (0) = -F n _2(0) and the initial conditions for n = and n = 1) 
we can compute the formal power series of the Fibonacci polynomial. 

> Fibonacci := proc(n,x) 

> if x=0 then sin(n*Pi/2) ~2 

> elif n=0 then 

> elif n=l then 1 

> else 'procname (args) ' 

> fi 

> end: 

> Fib := FormalPowerSeries (Fibonacci (n,x) , x) ; 

2 

Fib := sin(l/2 n Pi) (n - 1) (n + 1) 

/infinity \ 

' k (2 k) | 

\ (-1) (- 1/2 n + 1/2 + k) ! (1/2 n + 1/2 + k) ! x I 
) 1/ 

/ (n - 2 k - 1) (n + 2 k + 1) (2 k) ! I 
| 

\ k = / 
((- 1/2 n + 1/2)! (1/2 n + 1/2)!) 

/ infinity 

2 



n cos(l/2 n Pi) 



k (2 k + 1) 

\ (-1) (- 1/2 n + k) ! (1/2 n + k) ! x 
) 

/ (2 k + 1) ! 



\ k = / 

+ 1/2 

(- 1/2 n) ! (1/2 n) ! 

Note, that some of the factorials may have negative arguments, and hence the limits of the coeffi- 
cients must be considered. 

Let's for example compute Fiq(x) and Fn(x) which are polynomials of degree 9 and 10 respec- 
tively. For n = 10 = even we must only consider the second term of the solution Fib for k up 
to n/2 — 1 = 4 (as we shall show soon). Similarly for n = 11 = odd we only consider the even 
coefficients of Fib for k up to (n — l)/2 = 5. 

> limit (subs (Sum=sum, infinity=4, op(2,Fib)), n=10) ; 
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9 7 5 3 

x +8x +21x + 20 x +5x 

> limit (subs (Sum=sum, infinity=5, op(l,Fib)), n=ll) ; 

10 8 6 4 2 

x + 9 x + 28 x + 35 x + 15 x + 1 

In the second term of the above solution, the form of the odd coefficients are 
ncos 2 (fvr) (-l) fc (f + jfe)! (-§ + jfe)! 



(2 fc + i)! r(f-fe) 

where we used the identities 



2 (2* + l)! (§)! (-§)! 
ncos 2 (§vr) (-1)* r (l + § + k) T (l - (| - k)) 

2 (2fc + i)! r(i + f) r(i-f) 

ncos 2 (fvr) (-1)* r (1 + f + fc) T (f ) sin (f tt) 

2 (2fc + l)! r(l + f) r(f-/c)sin(^(f -fe)) 
cos 2 (fvr) (-i)fcr(l + §+fc) sin(fTr) 
(2fc + l)! r(f-fc) sin (7T (f — k)) 

cos 2 (| vr) r (1 + | + fc) = | (2fc +i)! C 4-t-V)i if n is even 

if n is odd 



r(x + 1) = xT(x) = x\ , 



r(x)T(i-x) = 



TT 



sin(-7ra;) 
and 

sin(a + b) = sin a cos b + cos a sin 6 . 
Similarly we get for the even coefficients 



sin 2 (§vr) (-l) fc (^(n+l) + fc)! n + 1 (-^n+l + fcy: „ 

(2k)\ n + 2/c + l (I( n + l))l n - 2fc - 1 (-±n + ±)l 

2 (n^) (_!)* (i(n - 1) + k)\ (-i(n + 1) + fc) ! 



sin 



(2fc)! (i(n-l))! (-|(n + l))! 

2 (n^ r (l + i(n - 1) + k) r (l - i(n + 



sin 



sin 



( 2k V- r(l + i(n-l)j r(l-i(n + l) 

2^) (_i)*r(l + i(n-l) + fc) r(i(n + l)) sin (f (n + 1)) 



( 2fc ) ! r(l + i(n-l)) r(i(n+l)-fc)sin((i(n + l)-A;)7r) 

2 (f7r) r(l + I(n-l) + fc ) = | ^If^l ifnisodd 

if n is even 



sin 



(2k)\ r (i( n + i)_^ 

Note that a 2 k+\ = for k > §, and a^k = for k > 



14 



If we make use of this additional information we end up with the following closed formula for 
the general Fibonacci polynomial F n (x). 



FJx) 



n/2— 1 

E o (2fc +i)i (^-Jiy x2k+l if n is even 

fe=0 (2fc)! (§("-!)-*)! It™ IS Odd 



Our implementation covers derivative rules and uses recurrence equations for the following 
families of special functions: the Fibonacci polynomials Fibonacci (n,x), the Bessel functions 
BesselJ(n,x) , BesselY(n,x) , Bessell (n,x) , and BesselK(n,x) (see |IJ, (9.1) and (9.6)), 
the Hankel functions Hankell (n,x) , and Hankel2(n,x) (see ]l|], (9.1)), the Kummer 
functions KummerM(a,b,x), and KummerU(a,b,x) (see Q, (13.4)), the Whittaker functions 
WhittakerM(n,m,x) , and WhittakerW(n,m,x) (see (13.4)), the associated Legendre functions 
LegendreP(a,b,x), and LegendreQ(a,b,x) (see Q, (8.5)), the orthogonal polynomials 
JacobiP(n, a,b ,x) , GegenbauerC(n,a,x) , ChebyshevT(n,x) , ChebyshevU(n,x) , LegendreP(n,x) , 
LaguerreL(n,a,x), and HermiteH(n,x) (see [U, (22.7) and (22.8)), and the iterated integrals of 
the complementary error function erfc(n,x) (see Q, (7.2)). 

As orthogonal polynomials are polynomials, for each fixed number n there is a simple DE of 
order one. This is the reason why we use different names from those in the packages orthopoly or 
combinat, as otherwise for fixed n evaluation occurs, and a first order DE is created which does 
not possess the structure of the polynomial system, and second, the derivative of orthopoly [T] 
can not be defined. For example 

> SimpleDE(LaguerreL(3,a,x) ,x) ; 

/ 2 \ 

I d I / d \ 

I F(x) I x + 3 F(x) + (a + 1 - x) I F(x) I = 

12 1 \ dx / 

\ dx / 

generates the second order DE which structurally characterizes the third Laguerre polynomial, 
whereas with the Laguerre polynomial out of the orthopoly package we get 

> SimpleDE (orthopoly [L] (3,a,x) ,x) ; 

2 2 
(18 + 15 a - 18 x + 3 a -6ax+3x) F(x) + 

2 2 3 2 2 3 

(6 + 11 a - 18 x + 6 a - 15 a x + 9 x +a -3a x+3ax -x) 

/ d \ 

I F(x) I = 

\ dx / 

We note that by a general result the sum, product, and composition with rational functions and 
rational powers of functions satisfying the type of DE considered, inherit this property fl^| , so that 
the algorithm is capable to generate the DE of a large class of functions. Here we give some more 
examples: 

> SimpleDE (exp(x/2)*WhittakerW (a, b, 1/x) ,x) ; 

43222 3 /d\ 

(x - 4x - 4x b + x + 4 a x - 1) F(x) - 4 x (x - 2) I F(x) I 

\ dx / 
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II \ 

Id 14 

+ 4 I F(x) I x = 

I 2 I 
\ dx / 

> SimpleDE(LegendreP(n,l/(l-x)) ,x) ; 
/ 2 \ 

2 I d I 3 / d \ 

x (x - 2) (x - 1) I F(x) I + n (n + 1) F(x) + 2 (x - 1) I F(x) I = 

12 1 \ dx / 

\ dx / 

> FormalPowerSeries(JacobiP(n,a,b,x) ,x=l) ; 

FPS/hypergeomRE: provided that -1 <= min(-l , -a-1) 

binomial (a + n, n) n a! 

/infinity \ 

' (- k) k| 

\ (- n + k) ! (b + a + n + k) ! (-2) (x - 1) I 
) 1/ 

/ (n - k) (a + k) ! k! I 
| 

\ k = / 
((- n) ! (b + a + n) !) 
Note, that here again for (— n + k)\/(—n)\ the corresponding limits must be considered. 



8 Asymptotic Series 

The same algorithm can also be used to compute asymptotic expansions. Note, that we only look 
for asymptotic series of the form of a Laurent-Puiseux series. Such series are unique, a property 
which is not given for general asymptotic expansions. As a consequence the results we compute 
may differ from the truncated asymptotic expansions returned by the Maple command asympt. 

One special thing of Laurent-Puiseux asymptotic expansions is, that they are only valid as long 
as the indeterminate approaches the expansion point from one side. If the expansion point is not oo, 
then one may specify with an option right or left from which side one approaches the expansion 
point. 

> FormalPowerSeries (erf (x) , x=inf inity) ; 



> FormalPowerSeries (arctan(l/x) , x=0, right); 



1/2 Pi 



/infinity \ 

' k (2 k + 1) 

\ (-1) x 
) 

/ 2 k + 1 



\ k = 

> FormalPowerSeries (exp (x) , x=-inf inity) ; 



/ 
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> FormalPowerSeries (exp(x) , x=inf inity) : 
FPS/FPS: ERROR: essential singularity 



> FormalPowerSeries (exp(x)*Ei(-x) + exp(-x)*Ei(x) , x=inf inity) ; 

/infinity \ 

\ (2 k + 2) 

) (1 + 2 k) ! (1/x) 

/ 



\ k = / 

> FormalPowerSeries (exp (x)* ( 1-erf (sqrt(x))) , x=inf inity) ; 
infinity 



(- k) (- k) (1/2 + k) 
\ (-1) (2k)! 4 (1/x) 
) 

/ k! 



k = 



1/2 



Pi 



By a plot of arctan(l/x) , e.g, one may realize that, indeed, the resulting series representation is 
one sided, only. 



9 Examples 

In this section we present some more results generated with the procedure FormalPowerSeries. 

First, we consider the FPS of the Airy wave function Ai. Since we use our own definition of the 
derivatives, we also must define the initial values for x = (see e.g. [PJ, (10.4.4)-(10.4.5)). 

> newAi(O) := 1/3" (2/3) /GAMMA (2/3) : 

> newAiPrime(O) := -1/3" (l/3)/GAMMA(l/3) : 

> FormalPowerSeries (newAi (x) ,x) ; 

/infinity \ 



1/3 



\ 



) 



(- k) (3 k) 
9 x 



■I 



1/3 



/ pochhammer(2/3, k) k! | 

\ k = / 
GAMMA (2/3) 

/ infinity 
\ 



1/6 

3 GAMMA (2/3) 



(- k) (3 k + 1) 
9 x 



) 



1/2 



/ pochhammer(4/3, k) k!| 

\ k = / 
Pi 

Note, that Maple is not yet able to compute even a truncated power series of Ai(x). The same 
holds obviously for the following example as long as the parameter a is left as an unknown. 
> FormalPowerSeries (x~a*sin(x~2) ,x) ; 
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infinity 

k (4 k + 2 + a) 

\ (-1) x 



) 



(1 + 2 k) ! 



/ 

k = 

The next two examples are interesting results which may be unexpected. It turns out, that 
both, e arccosh:r an d e arccosa; are f hyp er g eome t r i c type. 

> FormalPowerSeries (exp (arccos (x) ) , x) ; 

/ / k - 1 \ \ 



exp(l/2 Pi) 



infinity 

\ 

) 

/ 



\ j = 



(j + 1/4) 



k (2 k) 
4 x 



(2 k)! 



\ k = 



/ k - 1 



exp(l/2 Pi) 



infinity 
\ 

) 

/ 



I I 2 
i ! (1/2 + j + j ) 
I I 
I I 
\ j = 



k (2 k + 1) 
4 x 



/ 



(2 k + 1) ! 



\ k = 

> FormalPowerSeries (exp (arccosh(x)) ,x) ; 

/infinity 



- I 



\ 



) 



(- k) (2 k)| 

4 (2 k) ! x I 



\ k = 



(k!) (2 k - 1) 



+ x 



/ 



If we substitute Sum (the inert form of summation) by sum which tries to get a closed formula, then 



we get the well known identity e 



arccosh x 



X 



+ Vx 2 - 1: 



> eval(subs(Sum=sum, ")) ; 
In fact, both e arccosh:E an d e arccos:E are the special cases a = 1 and a 



2 1/2 
I (1 - x ) + x 



•1 of the function 



(x + Vx 2 - l) a . 

The following example has a logarithmic singularity at x = and hence the FPS of /' is 
computed. 

> FormalPowerSeries (arcsech(x) , x, real); 



ln(2) - ln(x) - 1/2 



/infinity \ 

' (- k) (2 k + 2) 

\ (1 + 2 k) ! 4 x 
) 

/ 2 

(k!) (k + 1) (2 k + 2) 

\ k = / 
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The next two examples are FPS expansions of two definite integrals, namely the polylogarithm 

X 

function and / evf(t)/tdt. Note that for the second integral we have used the inert function of Int 
o 

which is only a placeholder and does not try to compute a closed form. 
> FormalPowerSeries(dilog(l-x) , x) ; 

infinity 



\ x 
) 



(k + 1) 



(k + 1) 



/ 

k = 

> FormalPowerSeries(Int(erf (t)/t,t=0. .x) , x) ; 

infinity 



k (2 k + 1) 

\ (-1) x 
) 

/ 2 

k! (2 k + 1) 

k = 



Pi 



1/2 



The next two examples are FPS of functions which contain orthogonal polynomials. 

> FormalPowerSer ies (exp (-x) *LaguerreL (n , a , x) , x) ; 
FPS/hypergeomRE: provided that -1 <= min(-l , -a-1) 



binomial (n + a, n) a! 



/infinity \ 

' k k 

\ (-1) (n + a + k) ! x 
) 

/ (a + k) ! k! 



\ k = 
(n + a) ! 

> FormalPowerSeries(exp(-x~2)*HermiteH(n,x) ,x) ; 

/ infinity 



/ 



cos(l/2 n Pi) (n + 1) ! 



k (2 k) 

\ (1/2 n + 1/2 + k) ! (-4) x 
) 

/ (n + 1 + 2 k) (2 k) ! 



\ 



\ k = 
(1/2 n) ! (1/2 n + 1/2) ! 
/infinity 



/ 



sin(l/2 n Pi) (n + 1) ! 



k (2 k + 1) 

\ (-4) (1/2 n + k) ! x 
) 

/ (2 k + 1) ! 



\ 



/ 



\ k = 

+ 

(1/2 n + 1/2) ! (1/2 n) ! 

As we have seen above, we have three tools available to compute a formal power series repre- 
sentation of a given function, and it may happen for some examples, that several of these tools may 
be applied. Normally the solutions of a RE of hyper geometric type leads to the simplest results, 
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but this is not true in general. Hence we added the possibility for the user to choose a method by 
adding an additional argument which is either hypergeometric, explike or rational. 

> f := x*arctan(x)-l/2*ln(l+x~2) : 

> FormalPowerSeries (f , x, hypergeometric); 

/infinity \ 

k (2 k + 2) | 

\ (-1) x I 



1/2 



) 



/ (k + 1) (2 k + 1) I 
| 

\ k = / 

> FormalPowerSeries (f , x, rational); 

infinity 

k k (k + 2) 

\ ((- I) + I ) x 

) 1/2 

/ k k 

I (- I) (k +1) (k + 2) 

k = 

> f :=sin(x)*exp(x) : 

> FormalPowerSeries (f, x, explike); 



infinity 



\ 



) 



k 1/2 k 
(2 ) sin(l/4 k Pi) x 



/ 

k = 

> FormalPowerSeries (f, x, hypergeometric); 
/infinity 



k! 



k (- k) k (4 k + 1) 

\ (-1) 64 256 x 

) 



/ 



(4 k + 1) 



\ k = 



/ infinity 



k (- k) k (4 k + 2) 

! \ (-1) 64 256 x 

+ I ) 

I / 
| 

\ k = 



\ 



(2 k + 1) (4 k + 1) ! 



/ 



/ infinity 
| 

I \ 

+ I ) 

I / 
| 

\ k = 



k (- k) k (4 k + 3) I 
(-1) 64 256 x 



1 

(2 k + 1) (4 k + 1) ! (4 k + 3) I 

I 

/ 



10 Conclusion 



We have presented an algorithm, and its implementation in Maple, to compute FPS. Since it is 
a goal of Computer Algebra to work with formal objects, we think this is a very powerful and 
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valuable addition. 

The algorithm is beside of rational operations based on two basic tools: solving systems of 
equations and computing symbolic limits. No truncated series is ever computed. For the case of 
asymptotic series, one sided limits are computed and for all other cases complex ones. The power 
of the procedure for computing LPS stands and falls with the capabilities of the tool for computing 
limits. The algorithms which are used in Maple to compute limits are described in ||, [|. 

Further, in cases when the resulting RE cannot be solved explicitly, it can, in principle be used 
to calculate the coefficients iteratively in a lazy evaluation scheme. This is particularly efficient as 
the resulting RE always is homogeneous and linear, so that each coefficient can be calculated by 
finitely many of its predecessors. We note that the algorithm to find a simple DE works, and so such 
a RE always exists, if the input function is constructed by integration, differentiation, addition, 
multiplication, and the composition with rational functions or rational powers, from functions with 



the same property, see [|12[. This gives a huge class of functions to which the method can be applied. 
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